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Abstract 

We develop a discretized theory of thermal Casimir interactions to numerically calculate the interactions 
between fluctuating dielectrics. From a constrained partition function we derive a surface free energy, while 
handling divergences that depend on system size and discretization. We derive analytic results for parallel 
plate geometry in order to check the convergence of the numerical methods. We use the method to calculate 
vertical and lateral Casimir forces for a set of grooves. 
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I. INTRODUCTION 



Dispersion forces are long range interactions due to thermal or quantum fluctuation fields. The 
first theoretical derivation was due to Casimir [1] who predicted the attraction of two neutral plates. 
A general continuum theory, developed by Lifschitz, [2] addresses the non-additivity of these 
forces and is often used to interpret experiments [3]. Calculation of dispersion forces requires 
knowledge of the frequency dependent dielectric constant; information which is accessible from 
spectroscopy. The results of the theory can be expressed as a sum over Matsubara frequencies, 
[4]. However, one can isolate the zero Matsubara frequency and recognize that the corresponding 
contribution is both temperature dependent and independent of h; it is purely classical and de- 
pends on the static dielectric constant, e(uo = 0). The forces derived from all other frequencies 
depend on dynamic dipole fluctuations and require information on the frequency dependence of 
the dielectric permittivity. For biophysical systems, mainly composed of water and lipids, Ninham 
and Parsegian [5, 6] have shown that zero frequency gives a contribution to the interaction which 
is at least as important as the interactions coming from the UV region of e(uo). This feature makes 
biological materials rather unique, and justifies the study of the zero frequency contribution alone. 
Zero Matsubara frequency forces go under different names including static van der Waals forces 
[7], classical thermal Casimir forces [8], and Keesom interactions. For the simplest geometries 
one can calculate the free energy analytically [9, 10]. 

For more complex geometries this problem has been approached from various directions, both 
in its full quantum formulation, and in the high temperature, classical, regime. The range of 
techniques used is very wide, including Green's function formulations [11], world-line numerics 
[12], and path-integral formulations [13]. These techniques have allowed for great progress in 
linking the theory with experimental observation, but, as of now, no general method is available 
for arbitrary geometries and arbitrary values of the parameters. 

Recently, there has been a renewed interest in Casimir forces among soft condensed matter 
physicists as it was recognized that these forces can play an important role in biophysical sys- 
tems [7, 8]. Ninham and Parsegian observations have opened the road for new formulations of 
dispersion forces that focus exclusively on their classical part. Dean and Horgan showed how to 
calculate the free energy using a classical partition function without the machinery of Matsubara 
frequencies [8]. This formulation has the advantage of making the physics more transparent, as 
thermal effects are considered explicitly from the start, rather than obtained as a limiting case of a 
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quantum theory. 

In order to address a wider variety of geometries than is possible analytically, we introduce 
numerical methods that can be used to study thermal Casimir forces on a lattice. In particular 
we address the nature of the various terms appearing in the free energy of a dielectric system and 
study their variation with the discretization. This corresponds to the issue of regularization of 
divergences in a continuous infinite system. We apply our methods to a set of rectangular and 
sinusoidal grooves in order to characterize both longitudinal and transverse Casimir interactions. 

In this paper we are interested in computing the free energy due to thermal fluctuations of 
the electrostatic field between dielectric bodies. In this classical perspective, as noted above, the 
dielectric permittivity e is taken to be a function of space, e(r), but not a function of frequencies, i.e. 
e(r, u — 0). Our main goal is to develop a formalism which is a suitable starting point for efficient 
and versatile numerical methods, which can be applied to arbitrary geometries. Considering the 
classical, time-independent regime alone, allows us to present the general formulation and the 
numerical approaches on the simplest system. After a calibration which is run with convenient but 
non-physical values of e(r), we focus on systems with physically relevent material properties, that 
is systems for which the optical properties of the object immersed in water match those of water. 
This is the case if we choose a material whose refractive index is similar to that of water, n ~ 1.3, 
giving a dielectric constant for the material of t mat = 1.7, while €h 2 o = 80. 

In separate works we extend our approach to quantum regimes, considering the full dielectric 
function e(r, uo) [14]. We note also the existence of a second regime where thermal interactions 
can dominate even in the absence of matching of optical properties [15]; separations L satisfy 
L » he/ (2k bT). This case in only relevant when L » A/im and will not be discussed further 
here. 

The paper is organized as follows: In section II we discretize a set of dielectrics and derive the 
main theoretical results; in section III we study the discretized system using matrix diagonalization 
and Monte Carlo simulation. In section IV we present a continuum calculation which derives 
the free energy of a single homogeneous slab with periodic boundary conditions. This result is 
compared to the numerical results in section V, where we also show how the numerical methods 
can be used to calculate the free energy for grooved surfaces. 
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II. LATTICE FORMULATION 



The free energy of a dielectric in the absence of free charges is calculated from the partition 
function [16] 

Z= f V\D] n<HV-D)e-^, (1) 

J r 

with D(r) the electric displacement. The Boltzmann factor is the electrostatic energy in the pres- 
ence of a dielectric medium U = J (D 2 (r) /2e(r) ) d 3 r. Gauss' law V-D = is imposed in (1) as 
a constraint. In the continuum limit, and also in infinite volume, eq. (1) is ill-defined, as it contains 
divergences. We thus discretize in order to find an unambiguous definition of the free energy. We 
will then remove contributions to the free energy that diverge when the system size diverges, or 
mesh size is taken to zero, in order to calculate the remaining long ranged contributions to the free 
energy. 

We use a cubic lattice of iV nodes and 3iV links. We use the word lattice to refer to a regular 
grid. Since we are discretizing a macroscopic theory our lattice is not related to a physical atomic 
lattice. To be meaningful the lattice spacing should be large compared to atomic dimensions (so 
that one can discuss electrostatic interactions in terms of a dielectric constant) but small compared 
to the physical system at hand. 

For simplicity of notation we set the lattice spacing equal to I. Physical results depend on 
the ratio between the size of the bodies and their separation. Dimensional information can be 
reintroduced by standard scaling arguments. Physical quantities are associated with nodes or links; 
scalars "live" on nodes, vectors on links. Thus D nx is assigned to the link leaving the node n in 
the positive x direction. The discretized divergence V ■ D gives the total flux at a lattice site. The 



dielectric permittivity e is also associated with the links so that the discretized energy density takes 



£)2 £)2 ^2 

the form j^- + + j^. On a lattice, eq. (1) becomes 



2 f-f n 



links 



(2) 



Z = J V[B] ( n *(V-D) jexp 

\ nodes / 

where e h is the permittivity associated with the link I. Notice that, in periodic boundary conditions 
we impose N — 1 delta-functions; the constraint on the iV th node is automatically satisfied when 
N — \ are imposed since J d 3 r V-D = 

We impose Gauss' law by introducing an auxiliary scalar field = {4> n }n=i...N as a Lagrange 
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multiplier. It was shown in [17] that is the static electrostatic potential. Using the identity 



(2n) N Yl S(W ■ D) = J V[<f>\ exp 

nodes 



-i ^ v • D 



nodes 



(3) 



This identity involves the product over all N nodes, so we must introduce an extra constraint 5((f) N ) 
to remove the integral over <J) N . Dropping irrelevant prefactors. Eq. (2) becomes 



V[<j>]V\D] exp 



9 ^ ti 



links 



N 



<H<M H expH0V-D]. 



(4) 



nodes 



The scalar field <p is conjugate to V • D and is associated with the nodes. We integrate over the 
field D in eq. (4) and find an expression in terms of the field <\>. From here on we assume (3 — 1, 
or equivalently we absorb the temperature in the units of energy. 



\links 




V[<j>]5{(f> N ) exp 



links 



(5) 



Eq. (5) reduces the problem to that of calculating the partition function of a single scalar field with 
e-dependent gradient energy [8]. The only subtlety is the 5-function term in the measure. We will 
see below that its effect is to remove the integration over the constant mode, and to contribute an 
extra overall volume factor. 

Eq. (5) leads to a free energy containing several kinds of contribution. We wish to understand 
the lattice dependence of the various terms and remove those that diverge as we go to the contin- 
uum limit and scale as (V/s 3 ), where V is the volume and s the lattice spacing. What remains 
after taking care of these divergences corresponds to surface interactions [8, 10]. To separate these 
contributions we rescale to with x an arbitrary, positive function that, is associated with 
nodes. Then, from eq. (5) 

*=iK /2 n 



links 



nodes 



X 



1/2 



n d A 5 



\nodes 



>N 



exp 



links 



61 V 



Vx 



(6) 



We now focus our attention on the integral and on the remaining delta function. The first step is 
to recover simple Gaussian integrals by an orthogonal change of variable that makes the exponent 
diagonal. We write <f> n = a ni<^i, with n running over nodes, i = ... N — 1, and a ni are 
expansion coefficients. Then, the integral in eq. (6) becomes 




exp 



5 



fx* V 



CiiN 



(7) 
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where the A« are the eigenvalues of the operator (— -i=VeV^=) and are normalized eigenvectors. 
The lowest eigenvalue, A = 0, corresponds to the field configuration = N \/x, since for this 
configuration the gradient in the exponent acts on the constant vector. The normalized eigenvector 



is V 



/XWX2v 



/XJV 



TJT 



Thus, the A expansion coefficient of <p N is a 0N = ^/xn/ (X^ Xi) 



1/2 



We now integrate over a 

1 = ( e *) / (n ^) ex p [- A ^^] = ( e *) det * 

where the det* indicates the product over non-zero eigenvalues. Thus 

1/2 

E \ 1 (I<>1 



-1/2 



, (8) 




-1/2 



(9) 



\links / \nodes / \nodes 

We will now choose % to factorize the partition function into extensive terms and non extensive 
terms, separating the terms that diverge in the continuum limit from those that remain finite. 

We specialize to the case of interfaces between uniform dielectric materials. We discretize as 
follows. When a vertex belongs to a slab of material, all three links departing in the positive direc- 
tion from that vertex belong to the same slab. Thus, for a plane geometry, a slab of thickness a = 1 
contains one plane of L 2 vertices and 3L 2 links. This definition is not left-to-right symmetric, but 
it is convenient for periodic boundary conditions, as a slab of thickness a of e\ on a background of 
e is the same as a slab of thickness a' = L — a of e on a background of e\. We now chose x as 
Xn = \{tnx + tny + ^nz)- For a system where e is piece- wise constant this implies choosing x = e - 

From eq. (9) we find the free energy, T = — In Z 



+ - In det* 



i-v- £ vi- 



(10) 



T = -\j2 lnei + l E lne «-^ ln ( E ' 

links nodes \nodes 

The choice % = e, leads to a last term in eq. (10) which is homogeneous of degree zero in e. 
Thus scaling all values of the dielectric constant by a factor leave this contribution invariant. The 
other contributions contain e dependent terms which scale as N, or In N, giving divergences as 
V or In V in the infinite- volume limit. If we take the continuum limit keeping the volume finite, 
but sending the lattice spacing s to zero, these terms diverge as 1/s 3 or Ins. The determinant 
also has short-distance divergence in V, but this is independent of e and it is the same for all 
systems. It arises from the short-distance behavior of the operator — e -1 / 2 V • eVe -1 / 2 , which in 
all regions where e is a smooth function of r is the same as the short-distance behavior of the 
operator —V 2 . This contribution is the free energy of the vacuum. It can be subtracted defining 



Freg = F — \ m det* [—V 2 ] . In this paper we will be working with constant total volume, so this 
term can be ignored. We discuss these issues in detail in section IV. The surface free energy is the 
last term in eq. (10); it only contributes when e varies at an interface. We therefore define 



This interaction depends only on the ratio between the different dielectric constants of the uniform 
components of the physical system. For the water/lipid systems we wish to study this ratio is 
~ 1.7/80 ~ 50, we can than equivalently choose in our numerical codes e\ = 50 and e = 1. 
Expression (11) is free of short-distance divergences in any region where e(r) is smooth, since 
in those regions the short-distance behavior of the numerator and denominator of eq. (11) are the 
same. However, if e undergoes a sharp transition at a surface, there is an additional divergent con- 
tribution which can be ascribed to the self-energy of the interface between the two regions; a sharp 
interface affects modes of arbitrarily short wavelengths, therefore the large eigenvalue asymptotics 
of the numerator and denominator of eq. (11) do not cancel. This divergence is qualitatively dif- 
ferent from the volume divergences discussed earlier in this section. While the latter are the same 
in any continuum field theory and are usually dealt with using renormalization theory, the surface 
self-energy divergence is related to having defined the field theory on a singular background in 
which e(r) is discontinuous [18, 19]. In real physical systems treating the dielectric constant as 
discontinuous is an oversimplification; any transition between two material occurs over a small, 
but finite, distance. On the lattice, the lattice spacing gives a natural cut-off and the contributions 
coming from the surface self-energy scale as 1/s 2 . In this work we focus on the interaction en- 
ergy between separate surfaces, so we will subtract the surface self-energy by taking free energy 
differences between systems where the surfaces are rigidly translated but do not change shape. 

Eq. (10) is the fundamental equation we will use when extracting the surface interaction from 
the numerical results. We first determine the full partition function directly, either by matrix diag- 
onalization, or by Monte Carlo simulation. Using eq. (10) we extract the pure surface interaction 
eq. (1 1), by subtracting the extra terms in (10), which depend on the volume of the system and/or 
on the lattice. Explicitly 




(ID 




(12) 



where for simplicity we do not include the vacuum contribution from — V 2 . 
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III. NUMERICAL EVALUATION OF FREE ENERGIES 



We now compare two methods of evaluating the free energy, eq. (10), of the discretized system. 



The first method is a direct evaluation of the determinant det* 



-^V-eV^ 

Vx Vx 



-1/2 



The second 



is an implementation of the Monte Carlo algorithm, introduced in [20]. We consider a cubic box 
of side L, with periodic boundary conditions, with a cubic lattice of spacing s — 1. 



For small systems we evaluate the determinant det* 
methods. If we set x = U the partition function (9) then takes the form 



— Ly-eV^ 

Vx Vx 



-1/2 



Z= ([]f 1/2 ]i 3/2 def[-V-eV] 

\links / 



1/2 



with standard matrix 



(13) 



We write the exponent in (5) as a symmetric matrix acting on the field 0, J2 nodes ^M^^, where 
the non-zero elements of M are given by 

6 



nn=l 



Mi 



— e 



i,nn 



(14) 



where e^ nn indicates the value of the permittivity on the link connecting site % with the nearest 
neighbor in consideration. M is a symmetric matrix of dimension L 3 x L 3 . 

The determinant is calculated as the product of the eigenvalues of M, Aj. The free energy 
coming from det* [—V • eV] is 

1 i3 - 1 



where now the sum extends only from % — 1 to exclude the eigenvalue A = 0. The free energy, 
from (13) is then 



1 



det 



^lne-ln(L 3 ) 



Jinks 



(16) 



Comparing eq.(16) and(12), we see that the surface free energy is 



rf = Fdet - \ llle + \ ln ( E e ) - \HL 3 ). 

nodes \nodes / 



(17) 



The last thing we need to take care is the constant, e-independent contribution In det* (—V 2 ), 
which was discussed in section II, that represents the vacuum contribution (see also the Appendix 
for details). We can subtract this term by hand by computing it explicitly, or equivalently always 
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consider free energy differences. We adopt the second strategy, so that the constant factor InL 3 
entering in the last term of eq. (17) can also be ignored. 

As an alternative procedure, we have used the method of [20] to measure the free energy us- 
ing thermodynamic integration and Monte Carlo simulation. We sample the partition function (2) 
using a collective worm algorithm [21] to update the field D. We obtain free energy differences 
between two system: the one under study, and a reference system. In order to avoid surface diver- 
gences, the reference system has to be chosen as a rigid translation of the original system. In the 
example of parallel planar slabs, that we are going to analyze in detail, we take as "reference" the 
system where the surfaces are at maximum distance (i.e. L/2, in periodic boundary conditions). 

The free energy obtained by the simulation, T S i m , gives directly an evaluation of the left hand 
side of eq. (10). Using eq. (12), we conclude that the surface contribution is 

IV. ANALYTIC RESULTS 

The above discretization methods can be applied in arbitrary geometry. To calibrate them 
we will apply them to a flat slab, and compare the results with an analytic expression for the 
surface free energy, eq. (1 1), in the continuum limit. We consider a single-slab, piece- wise uniform 
system, periodic along the direction perpendicular to the slab, and infinite in the two transverse 
directions. Its geometry is shown in figure 1. Analytical results are well known for the similar 
system which is infinite in all three directions [8-10], but, to our knowledge, the result we derive 
here for the periodic system is new. 

Starting from the constrained partition function on the continuum, eq. (1), we derive the free 
energy of the periodic system following a procedure similar to [9, 10]. However, our derivation 
does not extract the classical result as a limit of a quantum system; it considers the classical system 
from the start. We believe that in such treatment the physics is more transparent and we are able 
to better control issues such as subtraction of divergences which, as we have seen, are important 
for the comparison of the analytical results with numerical data. In particular, since the final 
expression we derive is free of divergences, but the starting point, eq. (1), is not, we have to make 
sure that the finite quantities we calculate are consistent with the corresponding finite quantities 
defined on the lattice and computed through the numerical methods. We will show that expression 



5£M<) + i£M e )-5 ln (£< 

links nodes \nodes 



(18) 
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FIG. 1 : One slab geometry. 



(11) is essentially finite, (i.e. up to the -unphysical- divergent surface self-energy that appears in 
singular backgrounds, that will be discussed separately), so it is meaningful to compare it with 
numerical results. 

We will compute the free energy per unit transverse area. This quantity stays finite in the limit 
of infinite transverse size, and is defined, following eq. (11) 

1 det* r- e -V2 VeVe -i/2l 
^ = 2lf ln det*[-V3] (19) 



= _L {Tr* In [- e -V2 VeVe -i/2] _ Tr * ln [_ V 2] 1 ? (20) 
= T-T Q (21) 



where L t is the linear dimension of the system along the uniform directions, and the limit L t — > 00 
is understood. 

Determinants are cyclic invariant so that det(— e _1 VeV) is equivalent to 
det(— e~ 1 / 2 VeVe -1 / 2 ). In our test case, we find it more convenient to evaluate 

^=^|Tr*ln(-e- 1 VeV). (22) 

Let us consider the geometry in fig.(l), with e(r) piece-wise constant. We take our system to be 
periodic on [—§,§]. The surface of the slabs are perpendicular to the z axes. We separate the 
eigenvalue equation writing ip(x, y, z) = e lpxX+tpyV i/j z (z), with: 

e-\ep 2 -d z ed z )^ z = K^ z , (23) 

and p 2 = p 2 x + p 2 y . 

The eigenfunction along z is taken separately in the three regions 

^ a {z) = A a e^ + B a e-^ (24) 
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with a = 1, 2, 3. Integrating the eigenvalue equation over a small interval we find that ip satisfies 
the same boundary conditions as scalar electrostatic potential 0: 

1>.(z) = Mz), (25) 
e_dip-(z) = e + dip+(z), (26) 

which are derived by integrating (23) across the boundaries. Here, +, and — refer to left and right 
side of the interface. 

Eq.(23) gives p 2 + p 2 = A and p 2 + q 2 = A, from which we immediately deduce p — q. 
Moreover, we notice that A has to be real and positive. Indeed, consider an e(r) everywhere posi- 
tive, with the eigenvalues equation — (VeV)^ = eAip. Then, < J e|V^| 2 = — J ■?/>*( VeV)?/> — 
J V {ip*(eV4>)]. The last integral is zero because of the matching condition (26), leaving us with: 
< — / V ; *(VeV)V' = A J e\^\ 2 . As a consequence q has to be either real or purely imaginary, 
and in the latter case \q\ < p. 

Inserting the functions (24) in the eigenvalue equation (23) and using the boundary conditions, 
we find 

' e i - 6q \ 2 1 - cos(g(L - 2a)) _ ^ 
,ei + eo/ l-cos(gL) 
This equation determines the eigenvalues q n of the effective one-dimensional eigenvalue problem, 

eq. (23). The surface contribution to the free energy per unit area from (22) is[26]: 

T = ^Tr\n(-e-\eV t 2 + d z ed z )) 

= ^/ o d>p- 2 J2Mp 2 + V 2 n), (28) 

where q n 's are the solutions of eq. (27). The extra factor 2 in the integrand of the last expression 
comes from the consideration that the transformation q — > —q leaves both eq. (27), and the 
wavefunction (24), invariant up to a renaming of the coefficients. Summing over all solutions of 
(27) requires dividing by 2 to avoid double-counting. An equation similar to (28) holds for JF that 
we need to subtract, with the g n 's replaced by the appropriate eigenvalues for the uniform system. 

In the Appendix we show that the free energy, eq. (20), can be written in terms of a "spectral 
function" on the complex g-plane, Q(q): 

r~* = * ~ ^ = ^ [ *?l £ M/ + *»>fg, (29, 

where the integration contour 7 encloses all the complex solutions of the eigenvalue equation (27), 
for example the path shown in figure 2 (a), chosen in such a way that it avoids the branch cuts 
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(—zoo, —ip) and (ip, +ioo) from the logarithm. For the single-slab configuration, the appropriate 
function is: 

QW.!-^)' '-^-^ , (30) 
Vei + eo/ l-cos(gL) 

As shown in the Appendix, this choice corresponds to subtracting the energy of uniform empty 
space, the second term in eq. (20). 

To arrive at an explicit expression for eq. (29), we first rewrite the argument of the p-integral, 

/w ^* lV + a ||. (31) 

Rather than evaluating directly f(p), let us calculate its derivative: 

d , r dk 2p Q'(k) J dk 2p d. 

The integral along the large circle at infinity vanishes, and the only non-zero contribution comes 
from the paths along the imaginary axis. We deform the contour to two circles enclosing the two 
poles in ±ip, as shown in figure 2 (b). 

Eq. (32) is evaluated using the residue theorem: 

^ = ±lnQ(ip) + ^-\nQ(-ip), (33) 

whence 

f(p) = 2lnQ(ip) + 2f (34) 

where / is an integration constant, and we have used Q(q) = Q{—q)- Inserting eq. (34) in eq. (20), 
we obtain 

/*oo j: poo 

Fsurf = — / p\nQ(ip)dp+ — / pdp (35) 
4tt J 4tt J 

12 



The last term is still quadratically divergent. If we introduce a cut-off cu in momentum space, the 
last integral is O(fuo 2 ), the same scaling as a short-distance surface divergence. This is the explicit 
manifestation of the divergent surface self-energy described in section III, and it is caused by the 
singular background, characterized by a sharp transition in the dielectric profile. This divergence 
can be subtracted by setting the integration constant / = 0, as it does not depend on the dielectric 
constants. Notice however that the divergence will depend, in general, on the shape of the surfaces 
considered, so if we want to compare the free energies of systems whose interfaces have different 
geometries, we must treat this divergence more carefully. The self-energy can be interpreted as 
the finite tension of the surface interface, and one has to take it into account if one wishes to treat 
the surfaces dynamically. 

Setting the last term in eq. (35) to zero we find 

1 



Tsurf = ^ J p\nQ(ip)dp (36) 

Eq.(36) is gives us the free energy once we know the eigenvalue equation. In the appendix we 
show how the choices we made uniquely determine the function Q(q) as its zeros correspond to 
the solutions of the eigenvalue equation, and its poles correspond to the eigenvalues of the uniform 
system. These two requirements are equivalent to asking that the subtraction we performed to 
arrive at eq. (36) is exactly the second term in eq. (20). 
Using eq. (30), the surface free energy becomes 

2 



_ 1 f\m 1 ^i-^oV 1 — Ch[p(2a — L)] 



dp. (37) 



The integral over p is finite for separations a > 0; the integrand varies as pe~ 2ap for large p. Thus 
the surface free energy, eq. (11), is free of short-distance divergences and can be used in numerical 
computations. 



1 - " e 

ei+eo 



2 



dp. 



In the limit L — > oo, eq. (37) reduces to [8, 22] Fsurf = ^ pin 
The method of this section and the Appendix can be extended to systems composed of an arbitrary 
number of slabs through the use of transfer matrices [22]. 

V. NUMERICAL RESULTS 

We used both matrix diagonalization and Monte Carlo simulation to study the slab geometry 
of section IV in order to calibrate the numerical methods and evaluate discretization errors. We 
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FIG. 3: AFgurf as a function of surface separation from matrix diagonalization (solid) and Monte Carlo 
(points). The system size L = 12, dielectric constants eo = 1 and e\ = 2. The reference system is taken to 
be a system with maximal slab thickness, i.e. a? = L/2. 



firstly consider a volume V = N = 12 3 with lattice spacing s — 1. The dielectric constant of the 
slab is ei = 2 with background dielectric constant e = 1. These values of e do not correspond to 
any interesting physical system, and are far from the 1:50 ratio we seek when studying water/lipid 
systems. However one of the methods with which we wish to compare (Monte Carlo simulation) 
becomes inefficient for large dielectric contrasts. It is thus more convenient to work with reduced 
contrast in order to compare results. 

The matrix diagonalization is performed using the built in methods of Matlab. It requires about 
one minute on a 2GHz workstation to calculate the free energy as a function of the gap G, with G = 
1, . . . , 11. The Monte Carlo simulation was run with 40 points for the thermodynamic integration. 
For each gap measuring the free energy requires about one hour of simulation. Figure 3 shows that 
the agreement between the two methods is within error bars. 

The difference in computer resources needed in these two calculations of the free energy shows 
that Monte Carlo is not the best method for studying the large distance tails of the interaction 
between dielectric media. The surface interactions we want to measure are three to four orders of 
magnitude smaller than the self-energy ; their extraction requires extremely precise measures which 
are hard to achieve with statistical evaluations. We can estimate the scaling of the resources needed 
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in Monte Carlo as follows: The interaction between two interfaces scales as L 2 /G 2 = 0(1) which 
is to be compared to the volume free energy which scales as 0(V). Thus we require 0(V 2 ) Monte 
Carlo sweeps of cost 0(V) to generate statistically useful results for a single integration point. 
In addition, as the system size grows, more simulations points are needed for the thermodynamic 
integration. Considering all these factors, we see that the Monte Carlo method requires an effort 
that scales worse than L 9 . 

Matrix diagonalization is not affected by statistical noise. The complexity is dominated by the 
evaluation of the determinant: 0(V 3 ) = 0(L 9 ) with dense matrix routines. The real constraint for 
matrix diagonalization turns out to be the memory required for holding dense matrices. With 1GB 
memory available for computation, using standard routines, the largest system size we can consider 
is L max ~ 25; memory usage scales as 0(V 2 ) = 0(L 6 ). In order to study larger, more interesting 
systems we now specialize to objects which are translationally invariant in one direction. Because 
of the symmetry one can then use Fourier analysis to simplify the numerical problem. When we 
do this we find that we need only find the eigenvalues of a matrix of dimensions L 2 x L 2 instead 
of L 3 x L 3 . 

We write = ^ (f>^(f>^\ where (f) z 9 ^ are plane waves with q z = 2nn/L and n = 
0,1, ... ,L, corresponding to the eigenvalue 2(1 — cos(g 2 )) on a periodic lattice. Using this form 
for the eigenfunctions we find 

0M0 = ^0^M(g 2 )0^ (38) 
where the non-zero elements of the reduced matrix M(q z ) take the form 

4 

M i:i (q z ) = ^ e i.™« + 2e 2 (l - cos(g 2 )) 

nn=l 

M i:Tin (q z ) = -e i:Tin (39) 

and lndet* M = lndet* M(q z ). The sum over nearest neighbors runs only along the x and 
y axes, e 2 indicates the value of e along the uniform direction z. Using this approach, the new 
limit in size due to memory constraints becomes L max = Lflax ~ 125. However, in this case, 
computing time becomes the ultimate limiting factor because the determinant evaluation has to be 
performed L times, once for each q z , requiring an effort 0(L 7 ). A system size of about L = 50 
is the largest system we can consider on our same 2GHz processor, which leads to evaluating the 
determinant in about one hour. 
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gap 

FIG. 4: Ratio of Tth computed from the analytic formula eq. (37), and T m calculated from the matrix. The 
transverse dimensions are L x = L y = 60, while the longitudinal direction is L z = 30. 

We used this technique to compare the result obtained with matrix methods to our analytical 
calculations. Our analytic result eq.(37) is valid in the limit that the transverse directions of the 
system are large compared to the longitudinal direction. We thus study a system of dimensions 
60 x 60 x 30, where the smaller number refers to the longitudinal direction, and compare it to 
eq. (37) for e x = 50 and e = 1. Results are shown in figure 4 where we have plotted the ratio 
of the analytical result over numerical result as a function of surface separation. We see that the 
agreement is rather poor for the smallest separations, as might expected in a discretized system. 
For larger gaps the agreement significantly improves 

We now use the Fourier-matrix method to study the interaction between parallel grooves in di- 
electric surfaces. Using microfabrication techniques such grooves are easy to manufacture; recent 
theoretical studies have investigated the effects of corrugations on quantum Casimir free energies 
[11, 23]. We also note that one can study the Casimir energy as a function of lateral phase shift 
between the groove corrugations. Experimental evidence for lateral force was presented in [24]. 

For sinusoidal grooves analytic results were found using a second order perturbation in the am- 
plitude parameter while for rectangular grooves an exact numerical solution was found. The effect 
of corrugation wavelength and amplitude in relation to surface separation was investigated and it 
was found that for corrugation lengths much bigger than gap G, Proximity Force Approximation 
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(PFA) works well while in the limit of corrugation lengths much smaller than the amplitude a 
small distance regime is reached [25]. The first correction term to the parallel plate free energy, 
5 = T I? 'fi a t — 1, was determined. It was found that in the limit of the corrugation length A>G 
the correction 5 to the Casimir energy depends only on the corrugation amplitude, and it is pro- 
portional to In the opposite situation of A <C G the correction to the energy depends on both 
wavelength and amplitude and 5 ~ G~ 2 . 

Our numerical approach is valid for any groove profile, or indeed with any 2 + 1 dimensional 
interfacial profile, and for all separation regimes. We first consider the interaction of rectangular 
groove and a planar interface. We analyzed the dependence of the free energy on the corrugation 
wavelength A and on the surface separation G. In accordance to the analytical results for quantum 
Casimir, in the limit A > G we find the Proximity Force Approximation gives results in agree- 
ment with our numerical results. According to this widely used approximation, one considers the 
interaction of a small portion of the surface with a corresponding portion on the other surface 
which sits vertically above, assuming the surface is locally flat. The interaction between the two 
"plaquettes" is taken as the interaction of two parallel surfaces at the same vertical distance, and 
the contribution from all plaquettes is then added to obtain the total interaction. In a system of 
V = 40 3 , for the particular case of groove depth H = 4, and gap G — 6, at maximal wavelength 
A = 20 the agreement between the PFA prediction and the direct numerical result is around 90%. 
As A becomes smaller and comparable to the gap, the agreement drops to less than 20%. We have 
also investigated the free energy behavior in the limit A«G. We have indications that, as in the 
quantum system, the correction term 5 goes as l/G, but in order to be more certain one would 
need to study a larger system where ratios G/X ~ 10 2 can be reached. 

We now consider a sinusoidal groove, figure 5, with system size L = 46, amplitude of oscil- 
lation a = 6, wavelength L. The minimal distance G m i n between the grooves is reached when 
their relative phase is ir, and is equal to 1. In this configuration the vertical distance between the 
surfaces varies from G min = 1 to G max = 25. The shift S is measured in units of the lattice 
spacing as the displacement from the position where the grooves are perfectly aligned. Letting 
the top groove shift laterally over the bottom groove from the S = L/2 configuration (shown in 
figure 5) we measure the free energy as a function of the shift,jF(S'). J-'(S) gives information about 
lateral Casimir forces. When looking at the effect of the shift one expects the free energy to vary 
because the distribution of vertical distances changes as one groove is laterally displaced. Indeed 
when the shift is zero all the points of the opposing grooves have the same distance 2a + G min , 
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FIG. 5: Large sinusoidal grooves at S = L/2. The dashed profile indicates the position of the lower surface 
in the configuration S = 0. 

while for S = L/2, the distance between any two vertically aligned points varies from G min to 
Gmax = G min + 4a. But beside this projection of the vertical Casimir interaction along the direc- 
tion of the displacement, one expects to also detect lateral interactions between curved surfaces. 
Due to the non-additivity of Casimir forces, and because of the large deformation of the groove 
with respect to the flat geometry, these two contributions can not be disentangled. In figure 6 we 
present the results obtained with the reduced matrix diagonalization method, that gives the global 
interaction, and compare it to PFA, that only considers the effect of the local vertical displacement. 
PFA underestimates the free energy by factors that vary from around 10% when the grooves are 
mirroring (phase shift ir), up to 40% for phase shifts close to zero. This is what one would expect 
since PFA misses to consider collective effects, which are more important when larger positions 
of the surfaces are close. These violations of the proximity-force-approximation are not peculiar 
to the classical regime and similar effects are also found for the quantum system for significantly 
non-flat geometries [11, 12]. 
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FIG. 6: Free energy as a function of phase shift. 

VI. CONCLUSIONS 

We have presented a formulation of fluctuation induced interactions between dielectrics dis- 
cretized to a lattice. We analyzed how the contributions to the on-lattice partition function depend 
on the discretization and on the volume, and we have derived an appropriate definition for the sur- 
face interaction, which stays finite in the continuum limit. We compared two numerical methods 
to compute the free energy of surface interactions. The first is a direct evaluation of the matrix 
determinant. The second is a Monte Carlo simulation together with thermodynamic integration. 

While the constrained partition function allows the simulation of systems including full long- 
ranged Casimir interaction it turns out to be a rather inefficient method for extracting the asymp- 
totic interactions between bodies. The non-extensive surface interactions are rather easily lost in 
statistical noise. We consider that the use of matrix methods has considerable promise. Already 
interesting problems can be studied in the 2 + 1 dimensional in of translationly invariant systems, 
such as grooves, or blades. We note that we have only used the simplest dense matrix methods. We 
anticipate that the use of more specialized sparse matrix solvers should allow the study of larger 
physical systems, without the requirement of translation invariance in one dimension. 

Generalizations including the study of more elaborate dielectric functions, that include scale 
dependent dielectric response [16] as well as quantum effects will be presented in future papers. 
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APPENDIX A: REGULARIZATION OF ULTRA VIOLET DIVERGENCES 

In this Appendix we analyze the short-distance divergences in the parallel plate geometry dis- 
cussed in section IV. We address the subtraction of divergent vacuum energy and the divergent 
surface self-energy, showing that the latter arises from the singular nature of the background. 

In the continuum limit, there is always an e-independent divergent contribution in the free en- 
ergy which arises from the existence of modes of 4> with arbitrarily large-momentum. To see this, 
consider a system where e is everywhere uniform. The free energy corresponding to det*(— V 2 ) 
diverges as 



., , , „~n 2 . (A.l) 



n 



This divergence corresponds to large-momentum modes; it is present in the single-slab system 
analyzed in section IV, as the eigenvalues in eq. (27) also grow indefinitely. 
Consider now the free energy before the subtraction, eq. (28) 

^ = ^l2 Tr * ln (- e " 1VeV ) (A - 2) 

We now write the summation as a contour integral in the complex plane using Cauchy's theorem, 

where the spectral function Q(q) has simple zeros[27] for q = q n , and the integral is over an 
appropriate contour 7 that encloses all zeros of Q(q). We have used the fact that whenever Q(q) 
has a simple zero, Q'/Q has a simple pole with unit residue. Eq. (A.4) is clearly ill-defined, since 
the argument of the summation becomes arbitrarily large. 
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Notice that, if Q(q) has a pole of order p at q = q, eq. (A.4) is invalid, since close to q we have 

2^ „ (A . 5) 

(3(9) (9-9) 

and we get an extra contribution to the r.h.s. So, even formally, expression (A.4) holds only if 
the spectral function has zeros corresponding to the eigenvalues of the problem at hand, and no 
other zeros nor poles. However we can see from eq. (A. 5) that introducing simple poles q n in the 
spectral function is equivalent to subtracting the free energy T of a system that has eigenvalues 
q n . In this case therefore we formally have 

E ^p 2 + <£)-£ Hp 2 + Hp 2 + ^)f£y > ( A -6) 

where now Q(q) has q n as zeros and q m as poles. Depending on the properties of the spectral 
function, it can happen that the two terms on the l.h.s. are infinite, but the integral on the r.h.s. is 
finite, so that the expression (A.6), corresponding to the difference in free energies, is meaningful. 

Thus, the free energy difference between a given system and the vacuum can be written using 
eq. (A.6), in which 

Q(Q) = (A.7) 

Qo(q) 

where the zeros of Q(q) and of Q (q) are the solutions of the eigenvalue equations of the system 
under consideration and of the vacuum, respectively, and neither function has other zeros or poles. 
Using Eqs. (A. 3) and (A.6) we arrive at 

i.e. the equation (29) in section IV. 

For a uniform system in a periodic box of size L the eigenvalues are q n = ^p, so we can 
choose for example[28] 

Q (q) = cos(gL) - 1 (A.9) 

For the single-slab system considered in section IV, the eigenvalue equation (27) fixes Q(q) up to 
an overall constant, which can be in turn set by the requirement that Q(q) reduces to Qo(q) when 

Q(q) = ( (1 " cos[g(2a - L)]) - (1 - cos[gL]) . (A.10) 
\Cl + Co/ 

The resulting spectral function, 
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coincides with the one given in eq. (30), which as shown in section IV, when used in eq. (A. 8) 
gives a finite result for the surface free energy, up to a (sub-leading) divergence in the surface 
tension which will be discussed below. 

In section IV we found that we had to subtract an additional divergence that can be traced to 
the integration constant introduced in eq. (33). This originates from eq. (A. 8) being still ill- 
defined; only its derivative with respect to p is finite. In general however, there are situations 
when one can evaluate eq. (A. 8) directly and obtain a finite result. This is the case when when 
Q'(k)/Q(k) vanishes sufficiently fast for large \k\. Under this assumption, let us go back to the 
contour integration along the path shown in fig. 3 (a); we can neglect the integration over the large 
circle already at the level of eq. (A. 8), and integrate by parts in the remaining contribution along 
the cuts 

At this point, the branch cuts (—ioo, —ip) and (ip, —ioo) have disappeared, and we can again 
deform the contour to the two circles enclosing the two poles in ±ip, as shown in figure 2 (b). 
Then the residue theorem gives 

f(p)=2hiQ(ip), (A.13) 

which leads to the finite result, eq. (36). Notice that this procedure fails in the example of sec- 
tion IV, since Q(k) defined in eq. (30) does not have the required property for large \k\. 

The key observation is that, as we showed earlier in this appendix, Q(k) must be defined in 
such a way as to subtract the leading vacuum energy divergence, so that 

Q(k) = (A.14) 

where Qo(k) is a function giving the vacuum eigenvalue distribution, and Q(k) is such that firstly 
it gives the eigenvalue distribution of the system under consideration and secondly it reduces to 
Qo(k) for a uniform system. In all physical situations, the large momentum modes should always 
behave like in the vacuum beyond a certain threshold, determined by the properties of the material 
under investigation. Therefore it is natural to assume that in all physical systems, Q(k) — > 1 as 
| A; | — > oo, which automatically guarantees that Q'(k)/Q(k) —* for large \k\. In the case of a 
single plane slab, it is unphysical to assume that the dielectric constant changes discontinuously, 
because that would imply that all modes, with arbitrary short wavelength, are affected by the 
presence of the interface, as one can see from the matching conditions (26). It is more reasonable 
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to assume that modes with wavelengths shorter than a certain cut-off 5 (the molecular or atomic 
scale of the dielectric) will not be sensitive to the difference between the dielectric and the vacuum, 
therefore Q(k) ~ 1 for modes with \k\ > 1/5. Following these considerations it is clear that the 
divergence in the surface self-energy in section IV is exclusively related to the singular nature 
of the background under consideration, and that in a physical system with a smooth dielectric 
function e(r) the free energy in eq. (1 1) is finite. 
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